Microscopic theory of atom-molecule oscillations in a Bose-Einstein condensate 
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In a recent experiment at JILA [E.A. Donley et al., Nature (London) 417, 529 (2002)] an initially 
pure condensate of 85 Rb atoms was exposed to a specially designed time dependent magnetic field 
pulse in the vicinity of a Feshbach resonance. The production of new components of the gas as 
well as their oscillatory behavior have been reported. We apply a microscopic theory of the gas 
to identify these components and determine their physical properties. Our time dependent studies 
allow us to explain the observed dynamic evolution of all fractions, and to identify the physical 
relevance of the pulse shape. Based on ab initio predictions, our theory strongly supports the view 
that the experiments have produced a molecular condensate. 
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I. INTRODUCTION 

The subject of the coupling between atoms and molecules 
in Bose-Einstein condensates has attracted intense inter- 
est following recent experiments [1-4]. Particular exper- 
iments we shall focus on are those performed at JILA 
[3,4] where the strength of the inter-atomic potential of 
85 Rb was varied rapidly using specially designed mag- 
netic field pulses. This resulted in the loss of conden- 
sate atoms and the production of new components in the 
gas. One of these components is believed to be com- 
posed of molecules, and to be a molecular condensate. 
This is a remarkable achievement, with profound conse- 
quences for future work in the field. We shall show that 
this interpretation is fully supported by the theoretical 
work described in this article. We want to emphasize 
that the prediction of a molecular condensate arises nat- 
urally from the theory and does not have to be assumed 
at the outset. To make this prediction we use the micro- 
scopic theory of evolving condensed systems developed in 
Ref. [5] . This theory allows us to include the full dynam- 
ics of colliding pairs of atoms without the need for any 
assumptions about the nature of the states produced in 
the experiment. This theory gives us a generalization of 
the well known Gross-Pitaevskii equation (GPE) which 
includes the binary dynamics fully in the description of 
time evolving condensates. 

If the variation in the magnetic field occurs slowly in 
comparison with the duration of a collision one should 
expect to be able to use the standard Gross-Pitaevskii 
approach [6] to the problem. The derivation of the GPE, 
however, relies precisely on the assumption that colli- 
sions occur on a timescale small compared to all others 
in the problem [5]. This approximation, therefore, fails in 
this new experimental regime where the magnetic field, 
tuned in the region of a Feshbach resonance, varies on 
this timescale. 

As mentioned above the interpretation of the results 
of the experiment posits the production of bound molec- 
ular states, persisting at the end of the magnetic pulse 
sequence. Some theoretical treatments of the problem of 
molecules in condensates [7,8] separate out such states 



as a separate entity of the physical system at the outset 
of their calculation, i.e. physical observables associated 
with two-body bound states are described in terms of a 
molecular quantum field that emerges directly in a model 
Hamiltonian [9]. 

In this article we give a microscopic treatment of the 
evolution of a Bose-Einstein condensate in the presence 
of a time varying magnetic field that completely avoids 
assumptions on the nature of the states involved in the 
collision process. In fact we treat the binary events in- 
volving bound and free molecular states formed during 
the evolution in a unified manner. This is a most sensi- 
ble approach as strictly speaking, in the presence of the 
non-adiabatically time varying field there is no proper 
distinction between bound and free states. At the end of 
the pulse sequence we can of course resolve the final state 
of the gas into free and bound components. To do this 
we only use the assumption that the gas remains dilute 
and that binary encounters are the dominant collisional 
process. The evolution of pairs of particles from the con- 
densate into other free states or into bound molecular 
states comes from this treatment. Our theory strongly 
supports the view that the experiments have produced a 
molecular condensate. We should emphasize again that 
this conclusion comes from an ab initio prediction of the 
theory and not as an assumption. 

In the following sections we review briefly the micro- 
scopic theory of a dilute gas that we use in the analysis 
of the problem. We then show how the macroscopic evo- 
lution of the condensates is coupled in and out of the 
binary dynamics. We can then produce explicit expres- 
sions for the various components that are produced in 
the experiments. We have performed calculations, both 
for the case of the homogeneous gas and also for the case 
of a trapped condensate. The qualitative results of these 
two calculations agree but there are quantitative differ- 
ences that merit further study. Our results for the loss 
of condensate and the production of a heated component 
agree with those produced in the experiment. In addi- 
tion we are able to confirm the presence of a molecular 
condensate with all the physical properties we would ex- 
pect. 
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II. MICROSCOPIC DYNAMICS APPROACH 

A. Atomic mean field 

The microscopic dynamics approach [5] is based on the 
general many body Hamiltonian for identical bosons with 
a pair interaction V(r, t), 



H = J d 3 x^{x)H 1B { X )^(x) 



d 3 xd 3 y V t (x)^ t (y)^(x - y, *)V(y)V>(x). (1) 



Here, H 1B (x.) = -h 2 V 2 /2m + V tlap is the Hamilto- 
nian of a single atom containing the kinetic energy 
and the trapping potential. The field operators satisfy 
bosonic commutation relations, [^(xi), ip(x.2)] — 0, and 
[f/>(xi), ip^(x. 2 )] = <5(xi — x 2 ). In the situation studied 
in this article, the inter-atomic interaction is varied us- 
ing an external magnetic field pulse in the vicinity of a 
Feshbach resonance [4] and its time dependence is noted 
explicitly in Eq. (1). 

All physical properties of a gas of atoms can be de- 
termined from correlation functions, ( V^(x„) • • • V'(xi) )t, 
i.e. expectation values of normal ordered products of field 
operators with respect to the quantum state of the gas 
at time t. References [10,5] provide a general scheme to 
transform the exact infinite hierarchy of coupled dynamic 
equations for correlation functions into a more favorable 
form: The resulting equivalent set of dynamic equations 
for what are called non-commutative cumulants allows 
for a systematic truncation in accordance with Wick's 
theorem in statistical mechanics. In this article wc apply 
this truncation scheme to determine a closed set of equa- 
tions of motion for the relevant physical quantities. The 
derivation of the approach to the level of approximation 
required to study the phenomena reported in [4], i.e. the 
first order microscopic dynamics approach [5] , is given in 
Appendix A. 

The relevant physical quantities involve only the first 
and second order cumulants: 



*(x,t) = <V(x)) t , 

$(x, y ,t) - (V(y)V(x)) t - (V(y)MV>(x)} t , 

r(x,y,t) = (^(y)^(x)) ( - (^(y)) t <^(x)) t . 



(2) 



Here, ^(x, t) is the atomic mean field, $(x, y, t) the pair 
function, which plays an important role in the descrip- 
tion of correlated pairs of atoms, and r(x, y, t) is the 
one-body density matrix of the non-condensed fraction. 
The density of the gas at the position x and time t is thus 
given by n(x, t) = ( ^{x)ip{x) )t = T(x, x, t) + |*(x, t)\ 2 . 

In the first order microscopic dynamics approach the 
atomic mean field is determined through a closed nonlin- 
ear Schrodinger equation [5] : 



*ft^tf(x,t) = iTiB(x)tf(x,t) 

roc q 

- **(x,t) / dr* 2 (x,r)^-^(i,T) 



The collision term distinguishes the non-Markovian dy- 
namic Eq. (3) from the Gross-Pitaevskii approach and is 
determined through the coupling function 

h(t,r) = (2irh) 3 (O\V(t)U 2B (t,T)\O)0(t-T), (4) 

where U2 B {t,r) denotes the unitary time development 
operator of the relative motion of two atoms in free space, 
|0) is the zero momentum plane wave and 9{t — r) is 
the step function which gives unity for t > r and van- 
ishes elsewhere. Throughout this article the three di- 
mensional plane wave with momentum p is normalized 
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as ( r |p ) = exp(ip • r/K) / s/I-kTi . 



B. Non-condensed fraction 

In the first order microscopic dynamics approach the non- 
linear Schrodinger Eq. (3) determines not only the atomic 
mean field but also the pair function, and, in turn, 
the density matrix of the non-condensed fraction, V, in 
Eq. (2). The pair function is given by 

$(x,y,t) = -jf' dr J dVdVf(x',r)f(y',r) 



x-(x,y|^ p (t,r)|x',y'>, 



(5) 



where U 2 ® p (t, t) is the unitary time development oper- 
ator of two trapped atoms interacting through the pair 
potential V(t). The density matrix of the non-condensed 
fraction expressed in terms of the pair function is given 

by 



r(x,y,t) = J dV$(x,x',i)$*(y,x',i). (6) 

As shown in Appendix A, Eq. (6) assures both the posi- 
tivity of all occupation numbers and the conservation of 
the total number of atoms in the gas: 



/ 



d 3 x [|*(x,t)| 2 +r(x,x,t)] = N c (t) + N nc (t) = N. 

(7) 



The form of Eq. (6) suggests a separation of the num- 
ber of non-condensed atoms into a molecular fraction and 
correlated pairs of atoms after a time dependent magnetic 
field pulse of the kind reported in Ref. [4] as follows: The 
total number of non-condensed atoms is given by 

N nc (t) = J Ar(x,x,() = J d 3 xdV|$(x,x',t)| 2 

= / d 3 Rd 3 r\$(R,r,t)\ 2 , (8) 



to 



dr 



(3) 



where the position dependence was changed to two body 
center of mass and relative coordinates R = (x + y)/2 



2 



and r = y — x, respectively. Under the assumption that 
the trap is switched off at time ta n , immediately after the 
pulse, and the magnetic field is held constant at its final 
value, the energy states of the relative motion of a pair 
of atoms become stationary. A complete set of energy 
eigenstates is given through 



bv)(4>bv\ + 



d 3 P \4 +) )(^ +) \ 



(9) 



where 4>bv are the molecular bound states of the final 
pair potential and c/)^ are chosen as stationary scat- 
tering states which, at large relative distance, become a 
sum of an incoming plane wave with momentum p and 
an outgoing spherical wave (see, e.g., [11] or Appendix 
B). Replacing the spatial integration over the relative 
coordinate r in Eq. (8) in favor of the energy eigenstates 
in Eq. (9) the non-condensed fraction splits into a molec- 
ular part and a scattering part: 

N nc (t)= f d 3 R l(R,<M$W>l 2 

+ J d 3 P \(R,^Mt))\ 2 ], (io) 

where | $(t) ) = / d 3 Rd 3 r | R, r )$(R, r, t). The choice 
of eigenstates and the physical meaning of the contribu- 
tions to Eq. (10) depend on the experimental situation to 
be described. As will be shown in the next subsections 
the molecular part in Eq. (10) determines the number 
of atoms bound to molecules after the pulse while the 
scattering part describes pairs of atoms emitted from the 
condensate in a ballistic expansion. 



C. Molecular fraction 

The operator that determines the number of pairs of 
atoms in the specific bound state </>b in a gas with N 
atoms reads, in its first quantization form, 



1 N 



>b,ij 



(11) 



i,j = l 



where i and j indicate the pair of atoms. Expressed in 
terms of the atomic field operator Eq. (11) becomes: 

^ h = \ j d3Xl ^ X ' X ^ b<X ' 2 ~~ X l^b( X2 ~~ Xl ) 

x <5(±(x' 1 +x 2 )-±(x 1 +x 2 )) 

x V + (x'i )V f (x^)V(x 2 ) (12) 

The mean number of molecules in the state \4>b) is thus 
given by 

N h (t) = {M h ) t = \f d 3 r'd\d 3 RMr'W h (r) 
x (^(R+^ W t (R _^ (R r WR+ r } \ ) (13) 



where R and r are center of mass and relative co- 
ordinates, respectively. The fourth order correlation 
function in Eq. (13) can be factorized into cumulants 
(cf. Eq. (A4)), and truncated in accordance with the level 
of approximation of the first order microscopic dynamics 
approach: 



(V t (x 4 )V t (x3)V(x2)V(xi)) 4 

= <V t (x 4 )^ t (x3)) t <V'(x 2 )V(xi)) 4 



(14) 



The mean number of molecules in the state | <j>b ) can then 
be expressed in terms of a molecular mean field as 

JVb(t) = f d 3 R\* h (R,t)\ 2 , (15) 

where 

* b (R,t) = ^=J d 3 r<j>* h (r) 

x [$(R,r,t) + *(R+§,t)*(R- §,*)] . (16) 

The overlap of the molecular wave function <fib with the 
second, factorized term on the right hand side of Eq. (16) 
can be shown to be negligible in all applications described 
in this article. The molecular part on the right hand side 
of Eq. (10) is thus twice the number of dimer molecules 
in the gas, i.e. the number of atoms bound to dimer 
molecules. The wave function ^(R t), which yields the 
density of the molecular fraction, is thus obtained sys- 
tematically in terms of atomic field correlation functions. 
The derivation leading to Eqs. (15) and (16) does not de- 
pend on cf>b being a bound state. The number of pairs 
of atoms in any two body state is obtained in an analo- 
gous way. In Subsection II D we will apply an analogue 
of Eqs. (15) and (16) to determine the number of atoms 
emitted from the condensate during the magnetic field 
pulse. 



D. Burst of atoms 



In this subsection we show that the scattering part of the 
non condensed fraction on the right hand side of Eq. (10) 
determines the number of relatively hot atoms emitted 
in pairs from the condensate. To this end we consider a 
ballistic expansion of the gas at time ifi n , i.e. the trap is 
switched off and the magnetic field is held constant im- 
mediately at the end of a magnetic field pulse. This is 
illustrated schematically in Fig. 1. 

A relatively hot fraction, if present, will expand much 
faster than the condensate and can be detected far out- 
side the remnant condensate at a time t sufficiently long 
after ig n . For a gas with iV atoms the observable for 
the number of pairs of atoms with a relative coordinate 
between r and r + d 3 r, in its first quantization form, is 
given by J2i<j \ r ij ) ( r ij \d 3 r, where, as in Subsection II C, 
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i and j indicate the pair of atoms. In complete analogy to 
Subsection II C the mean number of pairs of atoms with 
a relative coordinate between r and r + d 3 r becomes 



n r (t)d 3 r = J d 3 R\^ r (R,,t)\ 2 d 3 r, 



(17) 



where 



* r (R,i) = -^[$(R,r,t)+*(R+|,t)*(R-§,t) 



(18) 




~P 




P 



FIG. 1. Scheme of a ballistic expansion after a magnetic 
field pulse. Pairs of atoms are emitted from the condensate 
with one-particle momenta p and —p. The center of mass 
of the pairs stays confined in the remnant condensate with 
a momentum spread determined by the spread of momenta 
in the condensate. At a sufficiently long time after the pulse 
a burst of relatively hot atoms can be detected outside the 
remnant condensate. 



p 2 /m. The spectrum of the pairs of comparatively hot 
atoms, i.e. the number of pairs of atoms with a relative 
energy in the interval E Te i . . . E Te i + dE le \ , is thus given 
by: 



n(E Te i)dE iei = -\ 



dn p / d 3 R\^ p (K) 



(21) 



where * P (R) = ( R, ^p +) |$(ifi n ) )/V% is the amplitude 
on the right hand side of Eq. (20) and dfl p denotes the 
angular component of d 3 p. The sum over all energy com- 
ponents of the spectrum yields 

n(E ml )dE Iel = ~ J d 3 P J d 3 R\(R, 0W|$(t fin ))| 2 . 

(22) 

A comparison between Eqs. (10), (15) and (22) shows 
that the total non-condensed fraction of the gas con- 
sists of molecules and a burst of comparatively hot atoms 
emitted in pairs from the condensate with a time of flight 
spectrum of relative energies given by Eq. (21). Whether 
the non condensed fraction becomes significant depends 
on the time dependence of the magnetic field, i.e. the way 
energy is released to the gas. 



III. DYNAMICS OF THE GAS 



At the relative distances r under consideration, which 
exceed by far the size of the remnant condensate, the 
second, factorized contribution to Eq. (18) is negligible. 
The energy spectrum of the relatively hot atoms can be 
obtained from an expansion of ty r in terms of the energy 
states in Eq. (9) that correspond to a release of the atoms 
from the trap: 

f r (R,<) = -^5>bv(rXR, 

+ 1 J d 3 p^(r)(R,cbMMt))- (19) 

The molecular wave functions in Eq. (19) have decayed 
at the relevant distances r that even exceed the extent of 
the remnant condensate. The corresponding molecular 
contribution to the right hand side of Eq. (19) is thus 
negligible. Taking into account that the scattering wave 
functions (j) p + ^ are energy eigenstates of a pair of atoms 
after the pulse the remaining part of the amplitude in 
Eq. (19) becomes, after a short calculation using Eq. (5), 



¥ r (R,f) = ^=J d 3 p^(r)(R, ^^(tfcO 
y( ,-i^(f-W/t 



(20) 



For two identical atoms the relative kinetic energy E re \ 
and the relative momentum p are related through E re \ = 



A. Feshbach resonance and magnetic field pulse 

In this Section we discuss the evolution of the gas when 
a specially designed homogeneous magnetic field pulse 
is applied to tune the inter-atomic interaction rapidly in 
the vicinity of a Feshbach resonance. Motivated by the 
experiment of Donley et al. [4] we study the time varia- 
tion of the magnetic field shown in Fig. 2. The magnetic 
field varies linearly in time within the subsequent time 
intervals. 




FIG. 2. Time dependence of the magnetic field. The field 
varies linearly within the subsequent time intervals. The nu- 
merical simulations in this article have been performed with 
the following pulse shape: Fall and rise times: 12 [is; hold 
time at B — 155.5 G: 13 /is; evolution time at B — 160 G: 
icvoive = 10...40/is; initial and final fields: B = 162.2 G. 
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A Feshbach resonance occurs when the energy of a 
bound state of a closed channel potential is tuned close to 
the dissociation threshold of the ground state potential 
[12]. This tuning of the interaction in the inter- atomic 
motion takes advantage of the Zeeman effect in the elec- 
tronic energy levels of the atoms. If the closed channel 
bound state approaches the threshold from below, the 
inter-atomic potential supports a shallow (metastable) s 
wave bound state. Around the resonance, a slight change 
in the energy difference of the potentials thus leads to a 
large variation of the scattering length. Neglecting the 
slow decay of the s wave bound state, the scattering 
length depends on the magnetic field through the rela- 
tion 

<B) = a hg (l - ^1^) , (23) 

where AB is the width of the Feshbach resonance and 
B is the resonant field. We consider the resonance of 
85 Rb at B = 154.9 G, with AB = 11.0 G [13], which has 
been used in [4]. For the background scattering length 
we use the value dbg = —450 aeohr [4], where aeohr is the 
Bohr-radius. 



performed for varying evolution times t voivo- The num- 
ber of atoms in each component showed an oscillatory 
dependence on £ C voivc with the frequency corresponding 
to the energy of the shallow two body s wave bound state 
in the evolution period. Moreover, a fraction of missing 
atoms was found oscillating at the same frequency. An 
interesting side result reported in [4] is that the visibility 
of the oscillations depended sensitively on the presence 
of the initial and final ramp very close to the resonance 
(at 155.5 G in Fig. 2). In the next subsections we shall 
explain these observations. 



B. Coupling function 

We will study first the coupling function of the non- 
Markovian non-linear Schrodingcr Eq. (3) for magnetic 
field variations as shown in Fig. 2. The coupling function 
h(t, r), given in Eq. (4), reflects the binary dynamics that 
enters the description of the condensate through Eq. (3). 
We will discuss to which extent the binary dynamics can 
already explain why the particular field pulse was needed 
to observe the oscillations between the condensate and 
the non condensed fraction of the gas. 




-4000 - X 

-8000 - \ 

150 155 160 165 170 

B [Gauss] 

FIG. 3. The scattering length a in units of the Bohr radius, 
as a function of the magnetic field B, in the vicinity of the 
Feshbach resonance. a(B) is determined from Eq. (23), using 
Bo = 154.9 G, AB = 11.0 G, and a bg = -450a Bohr . 

Figure 3 shows the scattering length as a function of 
the field B in the vicinity of the resonance for the values 
of B that are relevant in this article. At the initial time 
to, the magnetic field B = 162.2 G implies a scattering 
length of about 228ctBohr (cf. Fig. 3). The interactions 
then vary according to the pulse shape in Fig. 2. Simi- 
lar to the experimental procedure [4] we will determine 
the dynamic evolution of the gas for fixed time constants 
and field strengths of the initial and final pulses, but for 
different evolution times devolve • 

In the experiment [4] an adiabatic field variation fol- 
lowed the pulse sequence. Finally the trap and the mag- 
netic field were switched off, and the gas freely expanded 
before the number of atoms in the remnant condensate 
as well as a burst of relatively hot atoms were detected 
by absorption imaging. A series of measurements was 



15000 




FIG. 4. The coupling function h(t, r) as a function of t, 
describing the two-body dynamics driven by the magnetic 
field pulse in Fig. 2, for r = to, and an evolution time 

^evolve — 10 [AS. 



We have numerically determined the time dependence 
of h(t, t) in the two dimensional plane (t, r) for to < 
t < ttin and to < t < t using the methods described 
in Appendix B. Figure 4 shows h(t, t) as a function of 
t, for t = to and the time dependence of the magnetic 
field in Fig. 2, with t voivo = 10 /xs. In the evolution pe- 
riod between i 3 = 37 its and t 4 = 47 fis, h(t,r) oscillates 
with the frequency ^ cvo ivc — 200 kHz. This particular os- 
cillatory dependence is to be expected, as, according to 
Eq. (4), h(t,r) involves the two body time development 
operator, U%B{t, r): In this period the binary potential 
is stationary and supports a shallow s wave bound state 
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^evolve _ S p ec tral decomposition of Uz&(f, r) shows that 
the contribution of this bound state to h(t, r), within the 
evolution period, is given by: 



of the oscillation frequency of h(t, t) in the nonlinear dy- 
namics of the condensate described by Eq. (3). 



h(t,T) 



{2-Khf- 

x6(t- 



1 ^ (* ) I ^ V ° 1VC > < 0b VOlV ° I ^2B (* 3 , 

T)e~ iE -' 



°(t-t 3 )/h 



r)|0) 
(24) 



where ts is the initial time of the evolution period. The 
frequency of the oscillations in the coupling function thus 
corresponds to the bound state energy in the evolution 
period, i.e. ^evolve = \E^ volvc \/h. The amplitude and 
phase of these oscillations, however, depend on the time 
evolution before. 

The particularly large amplitude in Fig. 4 is achieved 
by the initial ramp close to the resonance at B = 155 G 
in Fig. 2. To illustrate the role of the first ramp, Fig. 5 
shows h(t,r) as a function of t, at r = to, for a trape- 
zoidal pulse. Here, the trapezoidal pulse is chosen sim- 
ilar to Fig. 2 except that the initial and final ramps to 
B = 155.5 G are cut off, i.e. the magnetic field is held 
constant at 160 G between t = 4 ^s and t = 80 [is. As 
shown in Figs. 4 and 5 the first ramp to B = 155.5 G 
causes a pronounced enhancement of the amplitude by 
a factor of about 20 for the optimized pulse in Fig. 2 as 
compared to the trapezoidal pulse. 
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FIG. 6. The coupling function h(t, r) as a function of t and 
r, for the same parameters as used in Fig. 4. 



FIG. 5. The coupling function h(t, r) as a function of t, for 
t — to, describing the two-body dynamics for a trapezoidal 
pulse, i.e. the magnetic field pulse in Fig. 2 but without the 
initial and final ramps to B — 155.5 G. The magnetic field is 
thus held constant at B = 160 G from t = 4 fjs to t = 80 fis. 



Figure 6 shows the real and imaginary part of h(t, r) in 
the two-dimensional plane (t, t), for the same parameters 
as in Fig. 4. The figure reveals that the amplitude of the 
oscillations during the evolution period between t = 37 yus 
and and t = 47 /is rapidly decays in r. The phase of the 
oscillations, however, is largely independent of r. A fur- 
ther analysis shows that these properties of the especially 
optimized pulse form in Fig. 2 assure the reappearance 



C. Homogeneous gas 

The dynamics of a homogeneous condensate driven by a 
magnetic field pulse of the form in Fig. 2 already exhibits 
all basic qualitative phenomena reported in [4]. We will 
therefore study the time evolution of the condensate as 
well as the final non-condensed fraction in detail for this 
idealized gas. Thereafter, we will discuss the corrections 
due to the presence of a trap in Subsection HID. All 
physical quantities under consideration are determined 
by the nonlinear Schrodingcr Eq. (3) , driven by coupling 
functions of the form of Fig. 6, with a variable evolution 
time t evo ive- 
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FIG. 7. The time dependence of the relative conden- 
sate fraction n c (t)/no remaining of the initial density of 
no — 3.9 x 10 12 cm -3 for three different evolution times 
tcvoive = 10 /us, 12 fjs, and 14 /is. The sequence of magnetic 
field pulses is chosen as in Fig. 2 with a field strength in the 
evolution period of f? evo ivc = 160 G. 



FIG. 8. The density of burst atoms (2xn(E re \) with n(E ie \) 
given by Eq. (21)) as a function of the relative energy E IC \ in a 
uniform gas for three different evolution times t cvo i V c = 11 /is, 
12 pis, and 13 pis. The external parameters are chosen as in 
Fig. 7. The atoms are emitted from the condensate in pairs 
with momenta p and —p. The relative energy E re \ = p 2 /m is 
related to the energy of a single atom, Eib = p 2 /2m, through 
Bib = Brei/2. The vertical line indicates the scale of the mean 
energies of burst atoms reported in Ref . [4] . 



Starting from a pure condensate with the density 
no = 3 x 10 12 cm -3 , Fig. 7 shows the relative atomic 
condensate density n c (t)/no, as a function of t, for three 
different evolution times i CV oivc = 10 fj&, 12 us, and 14 us. 
The initial conditions correspond roughly to the low den- 
sity measurements in Ref. [4] . After an initial loss period 
during to = < t < 35 us the condensate density shows 
a distinct oscillatory behavior around 80% of the initial 
density. The frequency of these oscillations very precisely 
matches the bound state frequency in the evolution pe- 
riod, i.e. z/evoive = \E^ olvc \/h = 200 kHz. After the 
evolution period, the second magnetic field pulse, which 
shifts the atoms in and out of the vicinity of the Feshbach 
resonance, causes the condensate fraction to develop to 
values between 55% and 85%. The final fraction depends 
on the phase of the intermediate oscillations at the end 
of the evolution period when the second resonant pulse 
starts. The remnant condensate density at time tg n , 
immediately after the pulse sequence, therefore, also os- 
cillates as a function of i 0V oivo- While the first ramp to 
B = 155.5 G in Fig. 2 drives the amplitude of the oscil- 
lations of n c (t) in t the second ramp in Fig. 2 amplifies 
the visibility of the oscillations in n c (ifi n ) as a function 
of tcvoive- The fast oscillations of the function h(t,r) in 
t at the very beginning and ending of the pulse sequence 
(see Figs. 4 and 6) have only a minor influence on the 
evolution of the condensate. 



The atomic mean field determines the pair function 
through Eq. (5) and, in turn, the molecular fraction in 
Eqs. (15) and (16) as well as the energy spectrum of com- 
paratively hot atoms, Eq. (21), after the pulse. Figure 8 
shows the density of atoms emitted in pairs from the con- 
densate as a function of their relative energy for a uniform 
gas under the conditions described in Fig. 7. As the mo- 
mentum spread of the center of mass motion of the pairs 
corresponds to the small spread of momenta in the atomic 
condensate the energy of a single atom in a pair is related 
to the relative energy through Eib = E re \/2. The spectra 
exhibit a damped oscillatory dependence on the energy 
with a first, dominant maximum below -EibAb = 150 
nK. The time of flight energy spectra in Fig. 8, as de- 
scribed in Subsection II D, do not correspond completely 
to the experimental procedure in [4]. The presence of a 
trap should also modify their shape in a noticeable way. 
We expect, however, that Fig. 8 reflects the typical en- 
ergy scales of the burst atoms in Ref. [4]. 

The experimental procedure did not allow for a direct 
detection of molecules. We thus identify the fraction of 
missing atoms reported in [4] as those atoms that are 
bound to dimcr molecules after the pulse sequence. The 
total density of unbound atoms is then given by the ini- 
tial density uq minus twice the density of dimcr molecules 
in the homogeneous gas. In the course of our studies we 
have determined the remnant condensate as well as the 
final molecular fraction as a function of the evolution 
time ^evolve from 10 to 40 us in steps of 1 us. The total 
length of the pulses has thus been varied between 84 /xs 
and 114 us. 

The results are summarized in Fig. 9. The solid line is 
an interpolation of the data for the remnant condensate 
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density relative to the initial density of no = 3.9 x 10 12 
cm~ 3 (filled circles) with the sinusoidal fit function pro- 
posed by Donley et al. [4]. The frequency of the oscil- 
lations corresponds to the binding energy |££ volvc |//i = 
|£ b (160G)|/ft = 200 kHz. In Fig. 9 the filled squares and 
their interpolation, i.e. the uppermost curve, show the 
fraction of atoms which are not bound to dimer molecules 
after the magnetic field pulse (Fig. 2). Number conser- 
vation allows to determine the density of the burst of 
relatively hot atoms directly from the total density of un- 
bound atoms and the remnant condensate (cf. Eqs. (7) 
and (10)). For this reason the dotted curve in Fig. 9, 
termed "burst of atoms" , has been obtained by subtract- 
ing the solid curve from the dashed curve. Both the frac- 
tion of "missing" atoms and the "burst" of atoms exhibit 
oscillations with the frequency ^ ovo ivc = \E^ volvc \/h in 
Evolve as well as the phase relation with respect to the 
remnant condensate reported in [4]. 



c" 0.4- 




study the influence of the trap and the inhomogeneous 
local densities of the gas on the relative magnitudes of 
these components. 



Gross Pitaevskii 
Gaussian 
Thomas Fermi 
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FIG. 10. The different initial states studied in a spherical 
harmonic trap with an oscillator length of about Iho = 3 ftm, 
i.e. a trap frequency of Uho — 80s -1 , corresponding to the 
geometric mean of the frequencies used in the experiment [4] 
(dJho = x/^radiai^axiaO- The atoms are exposed to a mag- 
netic field of B = 162.2 G which implies a scattering length 
of about 228 aeohr- Shown are the densities as functions of 
the radius corresponding to the exact ground state wave func- 
tion given by the time independent Gross-Pitaevskii equation 
(solid line), and a Gaussian (dashed line), which is the best 
least square fit to the Thomas-Fermi approximation (dotted 
line). All densities correspond to a number of 17100 con- 
densed atoms. 



FIG. 9. The remaining fraction of condensate atoms, 
n c (tfi n ), (solid line) together with the total density of unbound 
atoms (dashed line) in a homogeneous g function of the 

final time tfi n , elapsed after a magnetic field pulse of the form 
in Fig. 2. All densities are given relative to the initial density 
of no = 3.9 x 10 12 cm -3 . The external parameters are cho- 
sen as in Fig. 2. The filled circles and squares correspond to 
direct calculations of the remnant condensate and the molec- 
ular fraction. The solid and dashed curves are interpolations 
with the sinusoidal fit functions proposed in Ref. [4]. The 
dotted line indicates the "burst" of relatively hot (unbound) 
atoms emitted in pairs from the condensate as determined 
from the remnant condensate and the total fraction of un- 
bound atoms through number conservation (cf. Eqs. (7) and 
(10)). The filled diamonds correspond to direct calculations of 
the "burst" fraction obtained from integration of the spectra 
in Fig. 8 with respect to the energy. 



D. Trapped gas 

The studies of the homogeneous gas in Subsection III C 
allowed to identify the three components of the gas as 
observed in [4]. The purpose of this subsection is to 



We have performed these studies similar to Subsec- 
tion III C but in the presence of a trap. The trap has 
been idealized as a spherical harmonic oscillator poten- 
tial Vt r ap(r) = \mio\ Q r 1 with an oscillator frequency of 
about Who — 80 s -1 and a resulting oscillator length 
of about Iho — 3 /im. The trap parameters correspond 
to the geometric mean of the frequencies in Ref. [4], 

i.e. Who = \J t <-' 2 a diai 1 -''' axial- We have studied initially pure 
condensates of N = 17100 atoms with different local 
densities. The different initial density profiles are shown 
in Fig. 10. The solid line is the density obtained from the 
exact solution to the time independent Gross-Pitaevskii 
equation for 17100 atoms with a scattering length of 
about a(162.2G) = 228aj3ohr- This stationary conden- 
sate ground state exhibits a nearly perfect agreement 
with its Thomas Fermi approximation. The second state 
we have studied is a Gaussian (dashed line) which, for 
17100 atoms, is the best least square fit to the Thomas- 
Fermi approximation (dotted line). 
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FIG. 11. The remaining fraction of condensate atoms, 
A'c(tfin), (solid line), together with the molecular fraction 
(dashed line) and the burst of relatively hot unbound atoms 
(dotted line) , as a function of the total time ta n elapsed after 
the magnetic field pulse. The fractions are given relative to 
the initial number of atoms, i.e. iVo = 17100. The calcula- 
tions take into account the exact quantum dynamics of the 
trap for the initial state given, in Fig. 10, by the time inde- 
pendent Gross-Pitaevskii equation (a) and a Gaussian density 
distribution (b). 



Figure 11 shows the condensate, molecular and burst 
fractions at the end of the pulse (Fig. 2) for the exact 
solution of the time independent Gross-Pitaevskii equa- 
tion (a) and the Gaussian (b). The evolution times 
range between 10 [is and 17 /j,s such that the total time 
of the pulse ranges between 84 /is and 91 /is. The anal- 
ysis of the trapped gas does not show any qualitative 
differences from its analogue in a homogeneous gas in 
Fig. 9. All components oscillate with the same frequency 
v evolve = \E^ volvc \/h and exhibit phase shifts similar to 
those in Fig. 9. Figure 11 reveals, however, that the rel- 
ative magnitude of the components depends sensitively 
on the local densities in the initial inhomogeneous con- 
densate. 

To study the role of the trap potential and the one 
body kinetic energy we have performed the same analysis 
as in Fig. 11 in the local density approximation which ac- 
counts only for the dynamics of the nonlinear Schrodingcr 
Eq. (3) for a uniform gas but with different densities, 
weighted according to the initial condensate wave func- 
tion. The approximation, however, neglects the trapping 
potential and the one body kinetic energy. Figure 12 
shows the different fractions of the gas under the same 
conditions as in Fig. 11. A comparison of Figs. 11 and 
12 reveals no qualitative differences, but a quantitative 
dependence of all fractions of the gas on the trap po- 
tential and the one body kinetic energy. The oscillation 
frequencies are the same but the mean values and am- 
plitudes vary. In the local density approximation the 
fraction of burst atoms results considerably smaller than 
in the exact calculation. A further analysis of the time 



dependence of the condensate fraction, similar to Fig. 7, 
shows that the predominant influence of the trap poten- 
tial occurs during the initial ramp to B — 155.5 G in 
Fig. 2. In this period of the pulse sequence the scatter- 
ing length becomes comparable to the oscillator length 
[14]. This additional length scale is neither accounted 
for in the homogeneous gas nor in the local density ap- 
proximation. The results of this subsection show that a 
quantitative comparison of the magnitude of the differ- 
ent fractions in the gas with the experiment [4] should 
include the correct trap potential with the precise local 
densities in the initial condensate as well as the precise 
time dependence of the magnetic field pulse. 
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FIG. 12. The remaining fraction of condensate atoms, 
Af c (ifin), (solid line), together with the molecular fraction 
(dashed line) and burst of relatively hot unbound atoms (dot- 
ted line), as a function of the total time ifl n elapsed after 
the magnetic field pulse. The fractions are given relative to 
the initial number of atoms, i.e. TVo = 17100. The calcu- 
lations are performed in the local density approximation for 
the initial state given, in Fig. 10, by the time independent 
Gross-Pitaevskii equation (a) and a Gaussian density distri- 
bution (b). 



E. Interpretation of the results 

In the preceding sections we have analyzed the dynamics 
of an initially condensed Bose gas of 85 Rb atoms exposed 
to a magnetic field pulse of the form in Fig. 2. We have 
identified the different fractions of unbound atoms in the 
remnant condensate and a burst of comparatively 
hot atoms with kinetic energies of about Eyq/Ub > 100 
nK. We have further analyzed the dependence of the rel- 
evant physical observables on the precise external con- 
ditions. Our results strongly indicate that the fraction 
of missing atoms reported in [4] corresponds to atoms 
bound to dimcr molecules that could not be detected. 

Donley et al. [4] raised the question whether these 
molecules form a condensate. The analysis in Section 
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II provides a definite answer on the basis of the micro- 
scopic approach to the many body quantum dynamics: 
In Subsection II D we have predicted that the molecular 
fraction would stay confined in the atomic condensate in 
a ballistic expansion immediately after the pulse while 
the burst fraction of unbound pairs of atoms rapidly dis- 
perses. The expansion served as the first experimental 
technique to prove the presence of a condensate. 

The formation of a molecular condensate can be phys- 
ically understood from the nature of the external per- 
turbation of the initial atomic condensate: The mag- 
netic field pulse provides energy to form dimer molecules 
and correlated pairs of burst atoms but no momentum 
to drive their centers of mass. A further analysis of 
Eqs. (5), (16) and (20) shows that, indeed, the molecules 
as well as the correlated pairs of burst atoms exhibit 
the same momentum spread in their centers of mass as 
the condensate atoms. Momentum is transferred only to 
the relative coordinate of the burst atoms such that the 
total momentum of the pairs remains negligibly small. A 
typical density of a molecular condensate in a spherical 
trap is shown in Fig. 13. 



3e+ll r 



2.5e+ll - 



initial atomic density x 4/100 

molecular density 

final atomic density x 4/100 




r[nm] 

FIG. 13. The initial (dotted line) and final (dashed line) 
atomic condensate density, together with the molecular con- 
densate density, rib, (solid line) for tovoivc = 16 /is in a spheri- 
cal trap. The curves correspond to the calculations in Fig. 11 
(a). The atomic densities are multiplied by a factor of 0.04. 



As a great advantage of the microscopic dynamics ap- 
proach, the molecular condensate wave function, \I>b in 
Eq. (16), has been derived from the exact hierarchy of 
dynamic equations for correlation functions of field oper- 
ators with the general Hamiltonian Eq. (1) using a sys- 
tematic truncation scheme. Hence, the determination of 
the molecular condensate fraction does not rely upon any 
assumption on the existence of a molecular order param- 
eter. The analysis in this subsection is independent of 
the precise trap geometry or the density profile of the 
initial state. 
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APPENDIX A: THE MICROSCOPIC QUANTUM 
DYNAMICS APPROACH 



In this appendix we will derive the first order microscopic 
dynamics approach to Bose condensed gases. The present 
derivation includes the explicit dynamics of the non- 
condensed fraction which is not considered in Ref. [5]. 



1. Dynamic equations for cumulants 

In general, the quantum state of a gas at time t is de- 
scribed by a statistical operator p(t) and the correspond- 
ing expectation value of an operator O reads: 



(0) t = Tr[p(t)0] 



(Al) 



Quantum expectation values of normal ordered products 
of field operators are termed correlation functions. The 
exact dynamics of the infinite hierarchy of correlation 
functions is given by 



d 

i^(^(x B )--^(xi)) t 



([^ t (x„).-^(x 1 ),H]) t , 

(A2) 



where n denotes the number of field operators in the nor- 
mal ordered product and H is the quite general Hamil- 
tonian in Eq. (1). The set of Eqs. (A2) is equivalent to 
the many body Schrodinger equation and determines all 
physical properties of a gas of atoms. Due to the inter- 
action term in Eq. (1) any finite subsystem of Eqs. (A2) 
does not close, and any attempt of an approximate solu- 
tion relies upon a consistent way of truncation. A system- 
atic truncation scheme, based on cumulants of correlation 
functions, has been proposed in [10] for the dynamics of 
fermionic many body systems on short time scales. The 
extension of this work to interacting Bose gases [5] , with 
a modified truncation scheme, allows to describe the dy- 
namics on the time scales that are relevant in this article. 

The cumulants considered here are equivalent to the 
connected n-point functions and are usually defined as 
derivatives of a generating functional (cf., e.g., [15]): 



(V t (x„)-->(x 1 )) c 



<JJ(x„) <5J*( Xl ) 
x In ( el d; H/*(x)*/Hx)+J(x)</>t (x)] \ 



j=j*=o 



(A3) 
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The cumulants may also be derived recursively, and the 
first three orders read: 

<0i)H0i) c , 

(0102) = (010 2 ) C + (0l) C (02> C , 
(0 1 2 3 > = (0 1 2 3 ) C + (0l) C (0 2 ) C (0 3 ) C 

+ (0 1 ) C (0 2 3 ) C + (0 2 ) C (0 1 3 ) C 
+ (0 3 > C (010 2 ) C , 

; . (A4) 

For an ideal Bose gas, in the grand canonical thermal 
equilibrium, all cumulants containing more than two field 
operators vanish. This is a consequence of Wick's the- 
orem in statistical mechanics according to which every 
number conserving normal ordered correlation function 
can be expressed as a sum of products of all possible pair 
contractions conserving the operator ordering. Moreover, 
as the expectation value of a single field vanishes, the 
second order cumulants become (0i0 2 ) = (0i0 2 ) c - 
In accordance with Eq. (A4) the cumulants of an order 
higher than two vanish. In an interacting gas the higher 
order cumulants thus provide a measure for the deviation 
of the state of the gas from thermal equilibrium. 

For the applications discussed in this article, with 
a condensate present in the gas, the relevant cumu- 
lants contain; the non number conserving condensate 
wave function \I/(x, t) = (?/>(x))j, the pair function 
$(x, y, t) = (^)(y)^(x))j and the density matrix of the 
non condensed fraction r(x, y,t) — ( ^>t(y)^>(x) )£. The 
cumulant approach consists in transforming the exact hi- 
erarchy of dynamic Eqs. (A2) into an equivalent set of 
equations of motion for cumulants. The exact dynamic 
equations for the cumulants up to the second order are 
given explicitly in Ref. [5]. The truncation scheme of 
Ref. [5] consists in retaining, to the order of n, the exact 
dynamic equations for cumulants up to the order of n as 
well as the free time evolution of the cumulants of the 
order of n + 1 and n + 2. The free time evolution of the 
cumulants of the order of n + 1 and n + 2 is obtained 
by neglecting, in their dynamic equations, all products 
of cumulants containing n + 3 and n + 4 field operators. 

The first order microscopic dynamics approach (n = 1) 
results in a closed nonlinear Schrodinger equation for the 
mean field ^ that allows us to describe the dynamics of a 
condensate close to but also far away from thermal equi- 
librium [5]. The first order dynamic equations that lead 
to the non-linear Schrodinger equation read [5]: 

d 

*R^*(x,t) = ffi B (x)*(x,t) 

+ J d 3 yV(x-y,t)y*(y,t) 
x $(x,y,t) + *(x,t)*(y,t) 



Here the one and two-body Hamiltonians are denoted by 

Hib(x) - -ft 2 V 2 /2m + V4 rap (x), 
i? 2B (xi,x 2 ) = Fib(xi) +#ib(x 2 ) + V(xi -x 2 ,i). 

(A7) 

Equations (A5) and (A6) include the loss of condensate 
atoms into the non-condensed fraction as well as the back 
action of non condensed atoms on the condensate on time 
scales comparable to collisional durations. On longer 
time scales the non condensed fraction becomes dilute 
and its back action is neglected. 

The first order dynamics of the non-condensed frac- 
tion is determined through the conservation of the total 
number of atoms in the gas: 

f = | A[r(x,x,() + |f(x,()| 2 ] (A8) 

is a constant of motion. The corresponding approximate 
dynamic equation for T is then obtained from its exact 
counterpart in Ref. [5] through 

d f 

ih—T{K 1 ,x 2 ,t) = < i2i B (xi)r(xi,x 2 ,i) 

+ J d 3 yT/( Xl -y,t)<f>*(y,x 2 ,t) 



x $(xi,y,i) + *(xi,i)tf(y,t) 
- {xi <-> x 2 }*. 



(A9) 



ih— $(xi,x 2 ,i) = i? 2B (xi,x 2 )$(xi,X2,*) 
+ K(xi-x 2 ,t)tf(xi,t)tf(x 2 ,t). 



(A5) 
(A6) 



The density matrix of the non-condensed fraction 
r(xi,x 2 ,t), as given in Eq. (A9), is determined solely 
by the evolution of <J> from the initial time t up to the 
present time t. 



2. First order dynamics 

In this subsection wc will derive the non-linear 
Schrodinger Eq. (3) as well as Eqs. (5) and (6) for 
the pair function and the density matrix of the non- 
condensed fraction, respectively. For the purpose of solv- 
ing Eqs. (A6) and (A9) formally in terms of the mean 
field ^ it is convenient to change the representation from 
the configuration space to the one body energy states 
of the trap potential \<fii) or, for a homogeneous gas, 
into Fourier space. The corresponding single mode an- 
nihilation and creation operators obey the commutation 
relations [aj,aj] = Sij and the field operator becomes 
^( x ) — Si^i( x ) a j- I n tn i s new representation the cu- 
mulants up to the second order read: $i(t) — (aj)^ 
&ij(t) = (ajdi)t, and Tij(t) = (ajaj)£. We will abbre- 
viate the trap states \(j>i) by \i) in the following. 

In the new representation Eq. (A6) for the pair func- 
tion assumes the form 
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+ Yl < *' 3 I ^(*) | fei, k 2 ) [$ felfe2 (t) + * fcl (t)] , 
fci,fc 2 

(A10) 

where £"i is the eigenvalue of _Hib with respect to the 
mode function </>j(x). Equation (A10) can be solved for- 
mally in terms of the two body Green's function [16]: 

*«(*)= E [(^j|C/t 2 rap(^io)|fcl,fe)a> fclfc2 (to) 

fci ,k 2 



+ I dr(i,j\ G ( £ (t, t) V(t) I fei, )* fcl (r)* fe2 (r) 



to 



(All) 



Here is the retarded two-body Green's function, 



ihf t - H 2B (t) ) (t, r) = <J(t - r), ( A 1 2 ) 



which vanishes for t < t. The retarded Green's func- 
tion is related to the time development operator of two 
trapped interacting atoms, given by the time ordered ex- 
ponential 



through 



G^ ) (t,r) = l^-r)f/ t 2 4(i,r), 



(A13) 



(A14) 



where 9(t — r) is the step function that yields unity for 
t > t and vanishes elsewhere. In all applications in this 
article the gas is a dilute condensate at the initial time 
to- The initial pair function &k 1 k 2 (to) on the right hand 
side of Eq. (All) can then be neglected [5] and Eq. (All) 
thus yields Eq. (5) in position space. 

Inserting Eq. (All) into Eq. (A5) for the mean field 
leads to the closed non-linear Schrodinger equation 

+ E / dr(i,k 3 \T^(t,r)\h,k 2 ) 
fci,fc 2 ,fc 3 *° 

x * kl (T)* k2 (T)nM (A15) 

where denotes the retarded two-body transition ma- 
trix in the time domain: 



limit in the relative motion of two atoms can be per- 
formed in the collision term in Eq. (A15) [5]. The cou- 
pling function in Eq. (4) thus involves the time devel- 
opment operator of the relative motion of two atoms in 
free space, denoted by U 2 B(t, r). Transformed back to 
the position space Eq. (A15) then yields the nonlinear 
Schrodinger Eq. (3). 

In the first order microscopic dynamics approach the 
pair function, as given by Eq. (All) with 'f>h 1 k 2 (to) = 0, 
determines the density matrix of the non condensed frac- 
tion: Differentiation with respect to the time t shows that 



ry(t) = £$ ifc (t)$; fc (i) 



(A17) 



is the solution of Eq. (A9) for an initial pure condensate, 
i.e. Tij(to) = 0, which derives Eq. (6). The density ma- 
trix of the non-condensed fraction thus assumes the form 
of a partial trace over one coordinate of a two body pure 
state. The corresponding occupation numbers, i.e. the 
diagonal elements 



(A18) 



are positive, independent of the specific choice of the ba- 
sis set. 

As shown in Subsections II C and II D the non- 
condensed fraction consists of a molecular part and a 
burst of atoms emitted in pairs from the condensate 
with a comparatively fast relative motion. This separa- 
tion corresponds to a ballistic expansion of a gas that is 
released from a trap. The density of molecules in the 
bound state 4>b is described by the mean field \I>b in 
Eq. (16). In analogy to the collision term of the non- 
linear Schrodinger Eq. (A15) the molecular mean field 
can be expressed in terms of the atomic condensate wave 
function 

i r°° a 
ttb(R,*) = -^=jf dr* 2 (R,r) — h b (t,r). (A19) 

The corresponding coupling function, h\,{t,r), involves 
the overlap of the molecular bound state wave function 
4>h and the two body time development operator which, 
in all applications in this article, is excellently approxi- 
mated by the thermodynamic limit: 

h h (t,r) = {27Thf/ 2 (4> h \U 2B (t,T)\O)0(t-T). (A20) 

The energy spectrum of the burst atoms in Eq. (21) 
involves an amplitude ^ P (R), similar to Eq. (A19), ex- 
cept that the bound state <ph is replaced by the stationary 
scattering state </> p + ^ which is associated with the relative 
momentum p, i.e. 



T^\t,T) = V(t)S(t-T) + V(t)G { 2 + B ) (t,r)V(r). (A16) * p (R) 

As in all applications in this article the trap potential 
is slowly varying on the spatial scale determined by the 
range of the binary interaction V(t) the thermodynamic 



tfin 



dr I d 3 R'd 3 r' 



V2 J t 

x*(R' + r'/2, t)*(R' - r'/2, t) 
x| : (R,^ p +)|[/ t 2 4(t fin ,r)|R',r'). (A21) 
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Here, ifi n is the final time immediately after the pulse 
in Fig. 2. As <f>p(r) is not confined in space (see Ap- 
pendix B) the coupling function corresponding to v f p (R) 
should explicitly account for the discrete nature of the 
trap states also in the relative motion of two atoms. In 
Section III we have determined the energy spectrum of 
the burst atoms for a homogeneous gas, i.e. in the ab- 
sence of a trap potential. The coupling function of "fp 
then becomes similar to hb(t, r) in Eq. (A20). 



APPENDIX B: TWO-BODY DYNAMICS 

In Appendix A we have formulated the many body dy- 
namics of a condensed gas in terms of the unitary time 
evolution operator of two atoms interacting through their 
inter-atomic potential. In this appendix we provide a 
practical approach to determine the relevant low energy 
time evolution of two 85 Rb atoms that serves as an in- 
put to the microscopic dynamic description of a partially 
condensed gas exposed to the time dependent magnetic 
field discussed in Section III. The approach takes advan- 
tage of the fact that in all applications in this article the 
binary interaction of 85 Rb is dominated by the presence 
of a shallow s wave bound state. 



1. Resonance enhanced scattering 

In ultra-cold dilute gases the energies of two colliding 
atoms are usually sufficiently small for the differential 
cross sections to become isotropic. In accordance with 
effective range theory [16] the s wave scattering ampli- 
tude can then be expanded as 



f (k) = -a + ia(ka) + 0(k 2 ) 



1 + ika 



0(k 2 ), (Bl) 



where k is the wave number that is related to the rela- 
tive momentum of two colliding atoms through p = %k. 
Equation (Bl) effectively provides an expansion in terms 
of kl, where the first dominant length scale / is given by 
the s wave scattering length a. The next correction term 
involves the effective range of the binary potential V [16] , 
denoted as r c g in the following. In general, both a and 
r e ff depend sensitively on the detailed shape of V. 

When the binary potential supports a shallow s wave 
bound state the scattering length is positive and may 
by far exceed all the other length scales set by V. This 
situation is sometimes referred to as a zero energy res- 
onance [16]. The scattering amplitude is then given by 
— o/(l + ika), as obtained from the right hand side of 
Eq. (Bl), which corresponds to the contact potential [17]. 
Extending the collision energies p 2 jra into the complex 
plane Eq. (Bl) yields the T matrix [16] of the contact po- 
tential which assumes the separable, i.e. factorized, form: 



T 2 b(z) 



\x)Z(x\ 



1 + i^J mz/H 2 a 



(B2) 



where z = p 2 jra + ie is a complex energy variable and 
the complex square root is chosen with a positive imag- 
inary part. The wave function |x) and the amplitude £ 
are obtained as (r|%) = S(r) and £ = 4-Kh 2 a/m, respec- 
tively. The T matrix determines all eigenstates of the 
two body Hamiltonian. The pole on the right hand side 
of Eq. (B2) indicates that the contact potential with a 
positive scattering length supports a single s wave bound 
state with the binding energy 



E h = -?i 2 /ma 2 



(B3) 



The separable form of Eq. (B2) is quite general when- 
ever the T matrix is dominated by the pole of a shallow 
s wave bound state 0b • A spectral decomposition of the 
two body Hamiltonian then implies that at low collision 
energies the T matrix is well approximated by [16] 



\x)Z(x\ 



l-Z(x\G (z)\ X )- 



(B4) 



where 



\x) = v\M, 



(B5) 

and Gq(z) = (z + h A/to) -1 is the free energy dependent 
Green's function of the relative motion of two atoms. In 
few body scattering theory Eq. (B5) is usually referred to 
as the unitary pole approximation [18]. A further analy- 
sis shows that the T matrix in Eq. (B4) has a pole at the 
exact bound state energy of the potential V. The wave 
function x(r) accounts for the spatial extent of V . 

The long range behavior of inter-atomic potentials is 
determined by the van der Waals dispersion interac- 
tion Vvdw( r ) = —Cs/r 6 . The spatial extent of V is 
then characterized by the van der Waals length Z v dw = 
(mCe/h 2 ) 1 ^. As shown by Gribakin and Flambaum [19] 
the next correction to the binding energy for an inter- 
atomic potential modifies Eq. (B3) to [20]: 



E\, = —h 2 /m(a — a) z 



(B6) 



Here a is the mean scattering length given in terms of 
the T function through 



i r(3/4) 
a -27I vdW ^xW 



(B7) 



Figure 14 illustrates the dependence of the bound state 
energy of the shallow s wave bound state of the 85 Rb 
pair interaction on the magnetic field B. The binding 
energies obtained from Eq. (B6) arc sufficiently accurate 
to match a recent exact coupled channels scattering cal- 
culation [21]. Although the 85 Rb dimer, in the vicinity 
of the Feshbach resonance, is particularly weakly bound 
in comparison to usual molecular ground states Fig. 14 
exhibits a pronounced difference between Eqs. (B3) and 
(B6). This deviation from the zero energy resonance 
situation is related to the large van der Waals length of 
about l v dw = 164 aeohr- 
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FIG. 14. The energy of the uppermost s wave bound state 
below threshold of two S5 Rb atoms, Eh, as a function of the 
magnetic field strength B. The dotted line shows Eq. (B3) 
which corresponds to a zero energy resonance, with a(B) de- 
termined from Eq. (23), using B = 154.9 G, AB = 11.0 G, 
and cibg = —450 OBohr- The dashed line shows Eq. (B6) which 
accounts for the van der Waals interaction using Ce = 4660 
a.u. [13]. The dashed curve can be compared directly with 
a recent exact coupled channels scattering calculation [21]. 
The solid line is obtained from the binding energies of two 
parameter separable potentials with the rounded parameter 
rj = 5000 x ma Bohr /?i and the scattering length chosen in 
accordance with Eq. (23). 



2. The separable potential approach 

Direct insertion shows that Eq. (B4) exactly solves the 
Lippmann-Schwinger equation [16] 



T 2B (z) = V + VG (z)T 2B (z) 



(B8) 



as long as the actual potential V is replaced by the sep- 
arable potential 



Kep = |X>£<X|- 



(B9) 



Since the pioneering work of Lovelace [22] separable ex- 
pansions of potentials [23] have played an important role 
in nuclear few body physics, as they provide a systematic 
approach to solve two body scattering problems analyti- 
cally in a limited range of collision energies. In this sub- 
section we shall determine a separable potential of the 
form of Eq. (B9) that accurately describes the dynam- 
ics of two 85 Rb atoms in the relevant range of magnetic 
fields and collision energies. 

The unitary pole approximation in Eq. (B5) is obtained 
from spectral properties of the two body Hamiltonian 
and thus applies to inter-atomic potentials [24]. Equa- 
tion (B5) reproduces the exact bound state energy of V 
but the scattering length is only approximate. A fur- 
ther improvement can be achieved by choosing |x) and 
£ in such a way that the separable potential in Eq. (B9) 
matches both the energy of the last s wave bound state 
and the scattering length of V at the actual magnetic 



field. At the low collision momenta under consideration 
the corresponding plane wave states do not resolve the 
functional form of the wave function x(r). The specific 
form of x(r) is thus not relevant as long as xi r ) decays 
in space on the length scale set by the van der Waals 
length. We have chosen a Gaussian form which in mo- 
mentum space is given by 



X(p) = (PIX) = (2nh)-V 2 e-™ 2 ' 2mh . 



(B10) 



In position space \\) is of the form \{r) oc exp(— r 2 /2a 2 ) 
with a range parameter a = ^Jfvq/m. The separable 
potential is then parametrized by the two constants £ 
and r]. These parameters have to be determined at each 
magnetic field through matching the binding energy and 
scattering length of the separable potential to the val- 
ues of Eb and a of the actual interaction V. The bound 
state energy of the separable potential is the real energy 
z = E), at the pole of the T matrix in Eq. (B4), i.e. Eh 
is determined through 



l-Z(x\G {E b )\x)=0. 



(Bll) 



The scattering length is obtained from the zero energy 
limit of the T matrix as 



a = {2irhY 



;0|T 2B (0)|0) 



m (2^) 3 |(0| X )| 2 
47r?i 2 l/e-(x|Go(0)|x>' 



(B12) 



It turns out that the optimized parameter 77 is inde- 
penden t of B with the corresponding range parameter 
a = \Jfir\jm roughly given by ? v dw/2. In the applica- 
tions in this article we have used the rounded value of 
77 = 5000 x ma Bohr /h and determined £(B) in such a way 
that the separable potential matches exactly the depen- 
dence of the scattering length on the magnetic field in 
Eq. (23). The resulting dependence of the binding en- 
ergy of the separable potential on the magnetic field B is 
shown in Fig. 14. The separable potential approach, as 
proposed in this subsection, does not depend upon the 
accuracy of model potentials for all scattering channels 
and, instead, describes the complete low energy collision 
dynamics in terms of the scattering length a and the van 
der Waals dispersion coefficient Cq . Both parameters of 
the actual binary potential are accessible to experiment 
[25,13]. 

To illustrate the degree of accuracy of the proposed 
approach we shall provide a comparison of the low en- 
ergy scattering properties obtained from the separable 
potential with the exact solution of the Schrodinger equa- 
tion for a well known inter-atomic interaction. We have 
chosen the scattering of two ground state 4 He atoms for 
this comparison as the binary potential supports a single 
shallow s wave bound state and, within several decades 
of intensive study, all properties of the interaction have 
been determined very accurately, to a large extent, from 
first principle calculations [26] . The scattering length and 
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the binding energy have been determined recently from 
experiment in Ref. [27]. 

We have determined the separable potential, as given 
through Eqs. (B9) and (BIO), that corresponds to the 
Tang, Toennies and Yiu (TTY) 4 He interaction [26] with 
a = 188 ciBohr and C$ — 1.461 a.u.. For 4 He the ratio 
of the scattering length and the effective range, a/r c ff is 
of the order of magnitude of 14. Figure 15 shows the 
exact radial probability distribution of the bound state 
wave function 4 Hc 2 as well as the shallow s wave bound 
state of 85 Rb, corresponding to the separable potential 
approach at a magnetic field of 162.2 G, as obtained from 
the integral form of the Schrodinger equation 

|0b> = G (E h )V\<f> h ). (B13) 

The radius is given on a logarithmic scale. According 
to Eq. (B13) the asymptotic functional behavior of both 
bound state wave functions at far relative distances of 
the two atoms is determined solely through the free 
Green's function evaluated at the binding energy E\,. 
As their molecular states are very weakly bound both 
wave functions extend far outside the range of their pair 
interaction. 



of the inner well of the TTY potential in Fig. 15. 



— exact 




FIG. 16. The 4 He radial probability density corresponding 
to the s wave bound state of 4 He as obtained from the TTY 
potential [26] (solid line) and the separable potential approach 
(dashed line). 
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FIG. 15. The 4 He binary (TTY) potential [26], and the ra- 
dial probability densities corresponding to the shallow s wave 
bound states of 4 He as well as 85 Rb, as obtained from the 
separable potential, at a magnetic field strength of 162.2 G. 
The radius is given on a logarithmic scale. 



For 4 He the estimate corresponding to a zero energy 
resonance in Eq. (B3) gives \E^\/h — 25486 kHz while 
the formula of Gribakin and Flambaum in Eq. (B6) yields 
\Eb\/h — 26856 kHz. The exact binding energy of the 
TTY potential is \E b \/h = 27087 kHz. The comparison 
shows that E b and, in turn, the bound state wave func- 
tion, arc virtually completely determined by a and C§. 
Figure 16 compares the exact wave function of 4 He2 with 
the wave function obtained from the separable potential 
approach. The main small deviations occur in the region 



The T matrix determines the stationary scattering 
wave functions that correspond to the relative momen- 
tum p through the Lippmann-Schwinger equation [16] 



|<^ +) ) = \p) + G (p 2 /m + i0)T 2B (p 2 /m + i0)\p), 



(B14) 



where the energy arguments u p 2 /m + z0" indicate that 
the real energy p 2 /m is approached from the upper half 
of the complex plane. At low collision momenta p = hk 
the stationary scattering states assume the asymptotic 
form 



4+)(r) 



e tp - r/h + f (k) 



^ipr/h 



(B15) 



as soon as the relative distance r exceeds by far the 
range of the potential. Figure 17 compares the s wave 
scattering amplitude in Eq. (B15), for two 85 Rb atoms, 
obtained from the separable potential given by Eq. (B10), 
at B = 162.2 G, with the amplitude of the contact po- 
tential in Eq. (Bl). 
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FIG. 17. The real and imaginary part of the s-wave scat- 
tering amplitude fo(k) for 85 Rb. The solid lines are obtained 
with the two parameter separable potential, the dashed lines 
in the contact potential approximation, both with a scatter- 
ing length of a(162.2G) = 228aBohr- The wave number k is 
given on a logarithmic scale. 



FIG. 18. The real and imaginary part of the s-wave scat- 
tering amplitude fo(k) for 4 He. The solid lines are obtained 
with the two parameter separable potential, the dashed lines 
in the contact potential approximation. The dotted lines show 
the scattering amplitude of an improved contact potential ap- 
proach [28] that accounts for the effective range of the TTY 
potential. The bullets show the exact s wave scattering am- 
plitude for the TTY potential [26]. 



The pronounced deviations at k > 10~ 3 flg^ are 
related to the large van der Waals length of the 85 Rb 
interaction. This length scale is not accounted for by 
Eq. (Bl). The s wave scattering amplitude approaches 
Eq. (Bl) once the magnetic field is shifted further toward 
the Feshbach resonance at B — 154.9 G. The analogous 
comparison for 4 He in Fig. 18 may illustrate to which 
degree of accuracy the scattering from the long range 
part of the binary interaction is described by the separa- 
ble potential given by Eq. (BIO). In accordance with the 
small van der Waals length of helium of about 10 aeohr 
the deviations between the contact potential approach 
in Eq. (Bl) and the exact scattering amplitude are much 
less pronounced. Even the small deviations, however, 
are correctly accounted for in the separable potential 
approach up to wave numbers of about 3 x 10 _1 a^. 
The length scale related to this upper limit of the wave 
numbers roughly corresponds to the radius of the inner 
well of the TTY potential in Fig. 15. 



3. Dynamics 

In Subsection B 2 we have analyzed the static low en- 
ergy scattering properties of two 85 Rb atoms at a given 
magnetic field. In this subsection we shall determine the 
collision dynamics that enters the the many body theory 
of a partially condensed Bose gas through coupling func- 
tions of the form of Eq. (4). These coupling functions 
involve the complete unitary time evolution operator of 
two 85 Rb atoms, UzB(t,T), exposed to a magnetic field 
pulse as shown in Fig. 2. We shall apply the separable 
potential approach of Subsection B2 to determine the 
coupling functions as the effective low energy potential 
renders the time dependent Schrodinger equation into a 
practical form. 

We shall first determine the coupling function of the 
non-linear Schrodinger Eq. (3) denoted as h(t, r) in 
Eq. (4). The coupling function h(t,r) can be repre- 
sented in terms of the time developed zero momen- 
tum plane wave of the relative motion of two atoms, 
|C(t) > = U 2B (t,T)\0), in the form: 

h(t,r) = 9(t-T)(2irfr) 3 (0\V(t)\((t)). (B16) 

The wave function \((t)) is determined by the integral 
form of the time dependent Schrodinger equation through 



16 



IC(*)> 



01 



dT'G {t-r')V{T')\ar')), (B17) 



where Go(t) = 9(t)Uo(t)/iK is the two body Green's func- 
tion of the relative motion of two non- interacting atoms. 
To obtain the coupling function through Eq. (B16), 
on the basis of the actual binary potential V(t), the 
Schrodinger Eq. (B17) needs to be solved for all times 
(t, t) between the initial and final time of the magnetic 
field pulse and, moreover, at all relative distances in the 
argument of the wave function ({r,t). 
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FIG. 19. The variation of the s wave scattering length a 
corresponding to the magnetic field pulse in Fig. 2. 



The magnetic field pulse, however, releases a suffi- 
ciently small amount of energy to the gas that the actual 
potential V(t) in Eq. (B17) can be replaced by the effec- 
tive low energy potential V sep in Eq. (B9). Thereby, the 
amplitude £ = accounts for the time dependence of 
the magnetic field through the variation of the scatter- 
ing length a illustrated in Fig. 19. Equations (B16) and 
(B17) then yield the closed integral equation 



h(t,T) = (2nh) 3 \(0\ X )\ 2 at,T) 



+ W,t)J dr'(x\G (t-T')\x}h(r',r), 



(B18) 



where £(t,r) = £(t)6(t — r). The separable form of the 
effective low energy potential leads to a closed dynamic 
equation for h(t, r) which avoids to explicitly take into 
account the spatial dependence of C( r ^)- The coupling 
function h(t, r), as obtained from Eq. (B18), is shown in 
Figs. 4 and 6. 

The coupling function associated with the molecular 
condensate wave function \I/b is given in Appendix A by 
Eq. (A20) and denoted as h b (t,r). In Section III |* b | 2 
describes the density of 85 Rb 2 molecules at time ta n , im- 
mediately after the magnetic field pulse, in the bound 
state corresponding to the wave function in Fig. 15. The 
wave function £(r, tg n ) in Eq. (B17) determines /ib(ifim t) 
through: 

Mtfin.r) -^fin-r)(27r?i) 3 / 2 (^ b |C(tf i „)). (B19) 

Taking advantage of the separable form of the effec- 
tive low energy potential the Schrodinger Eq. (B17) in- 



serted into Eq. (B19) determines /ib(ifimT) in terms of 
the known coupling function h(t,r): 



tfln 



(2^) 3 / 2 (^ b |0) 



+ 



dt — <%> — M *' T) 



(B20) 



The molecular coupling function /ib(ifm, t), as a function 
of r, is shown in Fig. 20 for the magnetic field pulse in 
Fig. 2. The calculation of /ib(ifi n , t) has been performed 
with the 85 Rb2 wave function in Fig. 15 that corresponds 
to the shallow s wave bound state at the magnetic field 
of B = 162.2 G at the end of the pulse sequence. 
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FIG. 20. The molecular coupling function /i&(ifi n ,T), as a 
function of r, corresponding to the molecular wave function 
of two 85 Rb atoms at the final time of the pulse as depicted in 
Fig. 15. The magnetic field pulse corresponds to Fig. 2 with 

^evolve — 10 /^S. 



A relation, similar to Eq. (B20), with fa replaced by 
the stationary scattering state has been applied to 
calculate the spectral density of the pairs of burst atoms 
at a relative kinetic energy E Te \ = p 2 /m in a homoge- 
neous gas in Subsection III C. 
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